Arrhythmogenic influence of mutations in a myocyte-based computational model of the pulmonary vein sleeve

In the heart, electrophysiological dysregulation arises from defects at many biological levels (from point mutations in ion channel proteins to gross structural abnormalities). These defects disrupt the normal pattern of electrical activation, producing ectopic activity and reentrant arrhythmia. To interrogate mechanisms that link these primary biological defects to macroscopic electrophysiologic dysregulation most prior computational studies have utilized either (i) detailed models of myocyte ion channel dynamics at limited spatial scales, or (ii) homogenized models of action potential conduction that reproduce arrhythmic activity at tissue and organ levels. Here we apply our recent model (EMI), which integrates electrical activation and propagation across these scales, to study human atrial arrhythmias originating in the pulmonary vein (PV) sleeves. These small structures initiate most supraventricular arrhythmias and include pronounced myocyte-to-myocyte heterogeneities in ion channel expression and intercellular coupling. To test EMI’s cell-based architecture in this physiological context we asked whether ion channel mutations known to underlie atrial fibrillation are capable of initiating arrhythmogenic behavior via increased excitability or reentry in a schematic PV sleeve geometry. Our results illustrate that EMI’s improved spatial resolution can directly interrogate how electrophysiological changes at the individual myocyte level manifest in tissue and as arrhythmia in the PV sleeve.

www.nature.com/scientificreports/ model, but at the cost of a 15600-fold increase in the number of mesh blocks. Note that since convergence of the bidomain model is achieved at the mesh size ∼ x = 0.25 mm, further physiological details cannot be achieved by decreasing the mesh parameter using this model. Further comparisons of the bidomain model and the EMI model are found in 38 .
A range of genetic variants and mutations are known to increase the probability of developing AF, and a number of these occur in genes encoding the major cardiac ion channels [3][4][5] . In some cases, the electrophysiologic outcomes of these mutations have been well characterized, and provide a basis for asking whether a computational framework can be used to identify and characterize clinically meaningful arrhythmogenesis in the pulmonary vein sleeve. To implement such a framework using the EMI model, adequate data sets are available for the following mutations: N588K [39][40][41] , A545P 42 , E229V 43 , E375X 44 , A130V 45 , c.932delC 46 . Here, we investigate how each of these mutations change the likelihood of maintaining re-entry in a model of the specialized myocytes located at the transition from the pulmonary vein (PV) to the left atrial myocardium (LA). Furthermore, we examine how each mutation can alter the excitability of quiescent tissue in this region, and thus change the probability of triggering ectopic beats. In this numerical model, we introduce electrophysiologic heterogeneity in two ways: (1) by allowing each myocyte to be defined as a random combination of previously measured electrophysiologic phenotypes of PV and LA myocytes, and (2) by allowing electrical coupling to vary randomly among adjacent myocytes.
Our computations reveal that the effect of the selected mutations on the biomarkers of the action potential and cytosolic [ Ca 2+ ] varies significantly. None of the mutations markedly alter the resting membrane potential (RMP), whereas most have significant effects on the maximal velocity of the action potential (AP) upstroke, and AP duration (APD). As expected, we observed re-entry for a range of mutations that either slowed conduction velocity or shortened APD, namely N588K, E299V, A130V and c.932delC. Re-entry was not observed for A545P or E375X (and not for the baseline WT model). Tissue excitability was unchanged for most mutations, but significantly decreased for A130V, and significantly increased for c.932delC. Together, these results suggest that this high resolution (EMI) modeling approach can be applied to simulating the heterogeneous PV-LA junction, and thereby identify and discriminate the arrhythmogenic influence of clinically meaningful perturbations to myocyte electrophysiology.

Methods
In this section, we describe the models used in our simulations of the PV sleeve. First, we describe membrane models used to account for the ionic currents and intracellular ion concentrations of LA and PV cardiomyocytes. Next, we explain how each of the six considered mutations are represented in the models and the setup used for the finite element simulations of the EMI model. Finally, we provide the definitions of the biomarkers used in our computations.
Membrane models for PV and LA cardiomyocytes. To represent the resting membrane potential and action potential of both human LA and PV cardiomyocytes (CMs), we have utilized adjusted versions of the previously published base model 40,47 . Versions of our model have previously been used to simulate human induced pluripotent stem cell-derived CMs (hiPSC-CMs) and healthy adult human, canine, rabbit, guinea pig and zebrafish ventricular CMs 47,48 . In addition, this formalism has been used to study the short QT syndrome mutation N588K in ventricular CMs 40,49 . The base model used in this study is similar to the version used in 40,49 . Additional membrane currents (e.g., I Kur ) have been added to more accurately represent the AP waveform and underlying dynamics of atrial myocytes. The full base model formulation used here is found in the Supplementary Information. Unless otherwise specified, simulations characterizing the behavior of this model were performed at 1 Hz pacing frequency.
In order to represent the important differences between LA and PV myocytes, we have used published experimental data directly assessing differences in the membrane currents in myocytes isolated from both regions in the adult canine heart 22 . The only differences required to capture these effects in the LA and PV versions of the model are the maximum conductances of five transmembrane currents. Each has been shown to differ between LA and PV cells in 22 . The factors used to scale the conductance of the currents in the PV version of the model from the conductance of the currents in the LA version of the model are given in Table 1. The computational details of the human AP model used here are given in the Supplementary Information where the model also is compared to other AP models of the human atrial AP. have been linked to AF 3 . Five of these mutations alter the function of specific ion channels in the sarcolemmal membrane of CMs and one affects the function of gap junction elements (connexins) that couple neighbouring CMs. We assume the mutations affect the function of individual channel proteins, and that the function of these proteins and functional effects of the mutations do not differ between LA and PV myocytes. Therefore, the effect of the mutations is modeled in the exactly same manner for the LA and PV cases. Below, we list the selected mutations and explain how each is represented in the models. This information is summarized in Table 2. For all the cases, the models for the wild type (WT) version of the channels are unchanged.
N588K. The N588K gain-of-function mutation in the KCNH2 gene encoding the ion channels responsible for the rapidly activating delayed rectifier K + current, I Kr , is associated with short QT syndrome type 1 and other cardiac arrhythmias, including AF 41,50 . McPate at al. 39  A545P. The A545P gain-of-function mutation occurs in the KCND3 gene encoding the Kv4.3 alpha subunit of channels carrying the transient outward K + current I to , and is associated with AF 42 . A545P mutant channels exhibit increased peak current, and slowed current inactivation. In accordance with measurements published by Olesen et al. 42 , we represent the mutation by increasing the maximum I to conductance by 75% and the time constant of inactivation by 15%.
E299V. The E299V gain-of-function mutation in the KCNJ2 gene impacts ion channels responsible for the time-independent inwardly rectifying K + current, I K1 . This mutation is associated with short QT syndrome type 3 and increased incidence of arrhythmias, including AF 43 . To represent the mutation in the models, we use the WT and E299V versions of the atrial I K1 formulations provided in 43 . These formulations are based on the Grandi et al. human atrial myocyte model 52 and are fitted to WT and E299V I K1 measurements from 43 . The full formulation of the WT and E299V versions of the I K1 model are specified in the Supplementary Information.
E375X. The E375X mutation in the KCNA5 gene exerts loss-of-function effects in the Kv1.5 ion channels responsible for much of the measurable ultra-rapidly activating delayed rectifier K + current, I Kur , and has been causally linked to idiopathic lone AF 44 . Based on measurements of WT and E375X I Kur from 44 , we represent this mutation in our model by reducing the maximum conductance of I Kur to 10% of the WT value.
A130V. A130V is a loss-of-function mutation in SCN3B encoding the β 3 subunit of the Na + channel complex that is responsible for I Na , and was uniquely identified in a patient with lone AF from a study of 477 AF patients of Han Chinese descent 45 . Based on measurements of WT and A130V I Na from 45 , we represent this mutation in the model by reducing the maximum conductance of I Na to 25% of the WT value.
c.932delC. The c.932delC mutation in the GJA1 gene encoding the gap junction protein connexin 43 (Cx43) has been identified as a genetic mosaicism underlying lone AF in a patient for whom cardiomyocyte-specific expression of this variant was responsible for the pathology 46 . Based on measurement of the gap junction conductance for WT and the c.932delC mutation from 46 , we represent this mutation in the associated model by reducing the gap junction conductance ( G gap , see (2), below) to 25% of the WT value.
Finite element model of the PV 'sleeve'. We perform finite element simulations of a defined portion of the 'sleeve' of the PV using a model, referred to as the EMI model, which represents the geometry and electric potentials of the extracellular space (E), the cell membrane (M) and the intracellular space (I) (see, e.g., 28 www.nature.com/scientificreports/ for all myocytes k and neighboring myocytes k31 . Here, u e (in mV) is the electric potential in the extracellular space, e , and u k i (in mV) is the electric potential in the intracellular space of myocyte k, k i . Moreover, v k = u k i − u e (in mV) is the transmembrane potential of myocyte k, defined at the membrane Ŵ k = � e ∩ � k i , and w k = u k i − uk i (in mV) is the potential difference between myocyte k and its neighboring myocyte k , defined at the intercalated disc Ŵ k,k = � k i ∩ �k i . Furthermore, σ i and σ e (in mS/cm) is the conductivity of the extracellular and intracellular spaces, respectively, and C m and C g (in µF/cm 2 ) is the specific capacitance of the cell membrane and the intercalated discs, respectively. In addition, I k ion (in µA/cm 2 ) represents the sum of ionic currents through ion channels, pumps and exchangers expressed in the sarcolemma, represented by the base model described in the "Membrane models for PV and LA cardiomyocytes" section. The state variables of the base model are referred to as s k , and their dynamics are referred to as F k . Note that, I k gap (in µA/cm 2 ) is the current density between neighbouring myocytes, specified by the passive model where G gap (in mS/cm 2 ) is the conductance of the gap junctions and R gap (in k cm 2 ) is the corresponding resistance of the gap junctions. The parameters used in the EMI model simulations are specified in Table 3, and the EMI model equations are solved numerically using an MFEM 55,56 finite element implementation of the operator splitting algorithm introduced in 30,57 . A 500 ms simulation of cell to cell conduction within our PV sleeve substrate required a simulation time of approximately 8 days (190 hours).
PV sleeve geometry. We represent essential features of the geometry of the PV sleeve by constructing a collection of coupled myocytes that form a cylinder having a diameter of 1.5 cm (similar to the diameter of a human PV 58 ). Each myocyte is shaped as a cylinder with a diameter varying from 13 µ m at the cell ends to 14 µ m at the cell center, and each myocyte is about 120 µ m long 59 . In our formulation, the cylinder of myocytes consists of 393 cells placed longitudinally around the cylinder with 10 rows of cells positioned along the cylinder. The computational mesh of a single cell is illustrated in Fig. 1A and an associated cylinder of cells is illustrated in Fig. 1B Distribution of cell properties. We assume that the properties of the CMs forming the sleeve of the PV varies between those of PV and LA CMs but that there is no mean gradient from one end of the cylinder to the other. This arrangement is intended to specifically represent the subsection of the sleeve that is at the border between the vein and remote atrial myocardium. We represent this by drawing a number α that is 0 or 1 for each of the 393× 10 myocytes and letting the conductance of the five currents known to be different between PV and LA CMs (see Table 1) for each cell be given by where g LA and g PV are the LA and PV conductance values, respectively.
Somewhat similarly, the complex localized fibre arrangements and slow and complex conduction that have been observed in the PV sleeve 20 are represented by scaling the gap junction resistance between default values Table 3. Default parameter values used in the EMI model simulations, see, e.g., 30 . www.nature.com/scientificreports/ and values corresponding to a considerably reduced cell coupling. This is done by drawing random numbers β between 0 and 1 for each intercalated disc and letting where R * gap are the default gap junction resistances reported in Table 3. The numbers α and β are drawn randomly once from a uniform distribution and the same values are reused in all simulations.

Definition of excitability and re-entry.
We have performed EMI model simulations based on the PV sleeve to investigate how the selected ion channel or connexin mutations may affect the excitability and/or potential for re-entry in the tissue. Below, we describe the protocols used in these simulations.
Excitability. In order to investigate the excitability of the PV sleeve, we apply a 2 ms rectangular stimulus waveform ( I stim ) to 2 × 2 myocytes (see Fig. 1B). We evaluate how intense this depolarizing stimulus current must be to initiate an AP in the tissue, where a larger required I stim denotes reduced tissue excitability. We define excitability, E, as where I WT stim (in A/F) is the required stimulation current in the WT case and I MT stim (in A/F) is the required current for the mutation under consideration. Here, required refers to the smallest I stim capable of initiating an AP. www.nature.com/scientificreports/ Re-entry. In order to investigate how each selected mutation affects the likelihood of sustaining re-entry, we performed EMI model simulations in which 8 columns of myocytes around the cylinder are stimulated by setting their initial membrane potential to v = −10 mV. In addition, we enforce the condition ∂v ∂t = 0 in 4 columns of myocytes on one side of these 8 columns of myocytes, representing a refractory area. This condition is enforced until the propagating AP approaches these myocytes after traveling once around the cylinder. The approaching AP is monitored in each simulation, so that the point in time when the condition ends is adjusted for each individual simulation. We run the simulation for 500 ms and then determine whether the membrane potential of any of the myocytes is more depolarized than −65 mV anywhere in the cylinder at the end of the simulation. If so, we say that re-entry is obtained in the PV sleeve. The setup used to initiate re-entry is illustrated in Fig. 1C.
Defining biomarkers. In the results reported below, we utilize a number of different biomarkers, each computed from the solution of the ordinary differential equations defining the membrane models for LA and PV myocytes and from the EMI model simulations of the PV sleeve with a mix of LA and PV myocytes. In this section, we describe definitions of these biomarkers.

Electrophysiological. (a) Resting membrane potential (RMP)
The resting membrane potential (RMP) is defined as the membrane potential obtained 10 ms before the stimulation current that initiates the action potentials is applied.
(b) Action potential amplitude (APA) The action potential amplitude (APA) is defined as the difference between the RMP and the maximum value of the membrane potential obtained during an action potential.
(c) Maximal upstroke velocity ( dvdt max ) We define the maximal upstroke velocity ( dvdt max ) as the maximum value of the first derivative of the membrane potential with respect to time during the upstroke of the action potential.
(d) Action potential duration (APD) We define the action potential duration at p percent repolarization (e.g., APD50, APD80, APD90) as the time between the time when dvdt max is reached and the time during repolarization when membrane potential reaches a value below v min where v min is the minimum value of the membrane potential obtained between two action potentials and v max is the maximum value of the membrane potential during an action potential.
] transient amplitude (CaA) as the difference between the maximum value of the cytosolic Ca 2+ concentration obtained during an action potential and the lowest value obtained between two action potentials.
The Ca 2+ transient duration at p percent (e.g., CaD50, CaD90) is defined in the same manner as the APD.

Tissue dynamics. (a) Conduction velocity (CV)
We compute the conduction velocity in the EMI model simulations used to investigate re-entry by recording the times (after stimulation) at which the myocyte in row 5, column 50 ( t 1 ) and the myocyte in row 5, column 100 ( t 2 ) each first reach positive membrane potentials. The conduction velocity is then computed as where d = 50 · 120 µ m = 0.6 cm is the distance between the center of the two myocytes. This corresponds to the longitudinal conduction velocity of the cells, in the direction around the cylinder.
(b) Wavelength, For the excitation wave traveling around the cylinder shown in Figure 1, we also estimate the length around the cylinder that is covered by the wave. This conduction wavelength, , is computed as

Results
In this section, we present results from our simulations investigating the effects of selected ion channel protein mutations on the electrophysiological properties of a simulated functional sleeve of the PV. We first present key properties of the wild type models based mainly on the membrane currents of LA and PV myocytes, and we then examine how these properties are affected by each of the six mutations. Next, we investigate the effect of the mutations in EMI model simulations of the multicellular PV sleeve tissue. Specifically, we investigate whether the excitability of the CMs in the sleeve is altered and whether the mutations change the electrophysiological properties of the sleeve in a manner such that a re-entrant wave is sustained.
Properties of membrane models for LA and PV. In order to accurately represent the membrane currents and intracellular Ca 2+ fluxes of both LA and PV human myocytes and relate our findings to previous publications, we use a modified version of the base model from 40,47 . The differences between the LA and PV versions of our model are closely based on measured differences in transmembrane current densities from 22 . These are specified in Table 1. Figure 2 shows the AP and cytosolic [ Ca 2+ ] transient for the LA and PV versions of the base model. In Table 4, we compare a number of biomarkers for the two versions of the model. It is apparent that both human base models capture observed semi-quantitative differences between canine LA and PV myocytes from 22 . Specifically, the resting membrane potential is depolarized for the PV model compared to the LA model, and the maximum www.nature.com/scientificreports/ upstroke velocity of the action potential (AP) is considerably lower for PV than LA. In addition, the APD is longer for LA myocytes than for PV myocytes. All of these characteristics are consistent with the experimental observations 22 . Figure 3 shows the current densities for the different membrane currents of the LA and PV versions of the base model during an action potential. In Fig. 4, we also present the dynamic restitution properties of both models, which are key determinants of arrhythmia potential. The baseline difference in APD between the PV and LA models is retained across diastolic intervals (DI's), resulting in very similar restitution slopes, both of which are well below 1 until the DI approaches zero. Therefore, these restitution properties are unlikely to promote arrhythmogenic mechanisms involving APD alternans in either model. The effective refractory periods (ERPs) are more similar for the two models with the PV model still exhibiting a slight reduction. Again, the ERP restitution slopes are much less than unity until DI approaches only 10-20 ms (after 90% replarization) for the LA model.
Effect of the mutations on the membrane models. In Fig. 5 and Table 5, we investigate the effect of incorporating different mutations in the LA and PV versions of the base model as described in the "Modeling the AF-associated mutations" section. Note that a number of mutations (N588K, A545P and E299V) lead to a shortened APD for both PV and LA cells, and that the degree of shortening varies among the mutations. For the A130V mutation, we observe a significant decrease in the maximal upstroke velocity. For the E375X mutation, early phase II repolarization is markedly slowed, as is to be expected for this major loss-of-function I Kur mutant, although terminal repolarization and the indices of later APD are accelerated. Note that even though the E299V mutation affects the I K1 current, which is one of the most prominent currents during rest (see Fig. S3 in the Supplementary Information), the resting membrane potential of the LA and PV models appears to be relatively unaltered by this mutation (see Table 5). This may be explained by the fact that the difference between the WT and E299V I K1 is quite small for values of the membrane potential close to the resting membrane potential of the LA and PV models and more prominent for less depolarized values of the membrane potential (see Fig. S4 in the Supplementary Information). Since the c.932delC mutation affects only gap junction function, which is  Table 4. Biomarker values for the LA and PV versions of our base model at 1 Hz pacing. Selected biomarkers include the resting membrane potential (RMP), the action potential amplitude (APA), the maximal upstroke velocity of the action potential ( dvdt max ), the action potential durations at 90% and 50% repolarization (APD50 and APD90), the amplitude of the cytosolic Ca 2+ transient (CaA), and the cytosolic Ca 2+ transient durations at 90 and 50% (CaD50 and CaD90). These biomarkers are defined in the "Defining biomarkers" section.   . Action potential duration (APD) and effective refractory period (ERP) restitution curves for the LA and PV base models. For the APD restitution curve, the APD90 value is computed for a number of different diastolic interval (DI) values. The diastolic interval is defined as the time between the time point when the membrane potential has returned to 90% repolarization in the previous action potential (resulting from an S1 stimulation) to the time the next stimulation (S2) is applied. To compute the ERP restitution, a third stimulation (S3) is applied for each S2 DI. This S3 stimulation (lasting for 3 ms and of an amplitude about twice as large as was needed to generate an action potential from rest) is applied at progressively shorter time intervals after the S2 stimulation and the ERP is defined as the longest S2S3 interval that failed to capture 60 . Effect of the mutations on excitability. The effect of each selected mutation on the excitability of the myocytes in the PV sleeve has been investigated using the protocol described in the "Excitability" section. Figures 6 and 7 show the intracellular potential in these simulations at pre-determined, fixed points in time. In Fig. 6, we observe that a 2 ms stimulus current of 16 A/F is needed to initiate an AP in the WT tissue cylinder, and that this minimum stimulus intensity is unchanged for the N588K, A545P, E299V and E375X mutants. These findings are largely to be expected given that these mutations all impact K + currents, but did not influence resting membrane potential. In Fig. 7, we observe that for the c.932delC mutation, a smaller stimulus current is sufficient to generate a propagating AP, indicating an increased excitability for this mutation. On the other hand, a stronger stimulus is required for the A130V mutation, thus indicating decreased excitability.

RMP
Effect of the mutations on re-entry. To investigate how the same mutations affect the initiation and/ or maintenance of re-entry, we use the protocols described in the "Re-entry" section. Figures 8 and 9 show the intracellular potential of these simulations at some different points in time. In Figs. 8 and 9, we observe that in the WT case, a propagating wave begins to travel around the tissue cylinder, but upon reaching the stimulation site it encounters a refractory region and terminates. Thus, re-entry is prevented. The same mechanism is observed for the A545P and E375X mutations. www.nature.com/scientificreports/ In contrast, the shorter APs and refractory periods for the N588K and E299V mutations (Fig. 5), permit the approaching wavefront to re-excite this region. Consequently, the propagating wave continues to travel around the cylinder as a re-entrant wave.
For the A130V and c.932delC mutations, APD is not changed considerably (Fig. 5), but for A130V the maximal upstroke velocity is significantly decreased and for the c.932delC mutation, gap junction conductance is reduced by 75%. As a result, CV is considerably slower for both cases, thus permitting the myocytes at the stimulation site to fully repolarize, recover from refractoriness, and re-excite upon arrival of the propagating AP.
The wavelength of the conducted impulse. The wavelength of the conducted impulse is usually defined as the conduction velocity times the refractory period, see, e.g., 1 , and a re-entry wave is expected to be sustained if the geometric conduction pathway is larger than the wavelength. In our case, we expect re-entry to be possible if the wavelength is shorter that the circumference of the circle of myocytes illustrated in Fig. 1. In computations, the APD is often used as a proxy for the refractory period. In Table 6, we report the wavelength in our EMI model simulations, as defined by (7), in addition to the action potential duration (APD90), the conduction velocity (CV), the excitability (defined in the "Excitability" section) and the re-entry (as defined in the "Re-entry" section). In all but one case, we observe re-entry only if the wavelength is shorter than the circumference of the PV sleeve substrate. However, for the N588K mutation, we note that even though the wavelength is 4.9 cm and the circumference of the PV sleeve is 4.7 cm, re-entry is still observed. This is most likely because the APD90 value used in the calculation of the wavelength is only a rough approximation of the refractory period, and is discussed further in the Supplementary Information.

Summary of results.
To assess the ability of our EMI framework to distinguish mechanisms of arrhythmia in the PV sleeve we have imposed the effects of a range of arrhythmia-associated mutations on a schematic PV geometry incorporating heterogeneous intercellular coupling and cell electrophysiology at single myocyte resolution. The results of these simulations are, as mentioned above, summarized in Table 6. The simulations reveal that the dominant modes of arrhythmogenesis for these mutants are qualitatively consistent with their effects on single myocyte electrophysiology, but that they differ in tissue arrhythmogenic potential due to their severity or modulation of tissue electrotonic properties. To be specific, mutations that either considerably reduce APD (and refractory period, i.e. N588K and E299V), or slow conduction velocity (A130V, c.932delC) induce "leading circle" type reentry in our small and heterogeneous tissue geometry. Other mutations that cause less severe AP abbreviation, namely A545P and E375X, are unable to support genuine reentrant conduction in PV simulations, and thus are likely to require combination with additional AP shortening stimuli in vivo.
The two mutations that effect linear conduction velocity (A130V and c.932delC), also impacted excitability in the PV tissue by modulating its electrotonic properties. A130V reduced excitability by impairing the inward Na + current response to the trigger stimulus, while c.932delC enhanced excitability by reducing the electrotonic load imposed by the surrounding tissue. Table 5. Comparisons between computed WT and mutant biomarker values for simulations of the LA and PV versions of the membrane model driven at 1 Hz. These comparisons are based on changes in the same biomarkers as in Table 4, which are defined in the "Defining biomarkers" section. The numbers in parentheses report the deviation from WT. Note that the c.932delC mutation (in connexin 43) is not included here because these biomarkers are computed from simulations of the single myocyte LA and PV membrane models, and the c.932delC mutation does not affect any of the parameters of the membrane models.

Discussion
The borders between the PV sleeves and working atrial myocardium are small regions of tissue exhibiting pronounced electrophysiologic heterogeneity and structural complexity. Here we have asked whether a recently developed mathematical model (EMI), which is specifically built to interrogate heterogeneities at the level of single cardiac myocytes, is capable of identifying clinically meaningful arrhythmogenic tendencies in this region. We specifically selected a range of genetic mutations known to contribute significantly to the genetic basis of either lone AF, or AF with only comorbid arrhythmias (e.g., short QT syndrome). We then assessed whether the arrhythmogenic impacts of these mutants could be identified by an EMI implementation of a thin crosssection of the PV-LA border. Our model incorporates heterogeneity of myocyte electrophysiology as defined by patch-clamp measurements in isolated PV and LA myocytes, and sharp microelectrode tissue recordings. This approach was readily capable of assessing the proarrhythmic influence of these mutations, and also of revealing one dominant mode of arrhythmogenesis for each. Specifically, it was able to identify and begin to characterize the propensity for re-entrant behavior among mutants that expand the region of excitable tissue (in this small geometry), either by reducing conduction velocity or abbreviating the tissue APD (and hence the tissue refractory period). This computational approach was also able to distinguish the propensity for two mutations impacting I Na (A130V) and Cx43 (c.932delC) to markedly alter tissue excitability; and as expected, to do so in Figure 6. Excitability analysis for WT, and for the N588K, A545P, and E299V mutations. Intracellular potential of myocytes in the PV sleeve is illustrated at four different times (rows) for each mutation (columns). In all cases a regenerative AP was generated and propagated around the tissue cylinder. The stimulus intensity (applied to 2 × 2 cells for 2 ms) required to initiate the propagating action potential in each case is provided in the title for each mutation. Note that to improve the visibility, the transversal axis of the myocytes has been stretched by a factor of 20. The role of high-resolution modeling for understanding arrhythmogenesis in the PV sleeves. Since the initial identification of the PV-LA junctions as the primary sites of AF initiation 18 , much progress has been made towards understanding the unique characteristics of this region that so strongly predispose it to generating abberant and disruptive electrical activity. Some of these factors are structural, and include the geometry of the tissue (a long, thin cylinder of myocardium directly bordered by the insulating fibrous veinous wall), sharp discontinuities in myocyte fiber direction at the PV-LA border 61,62 , and characteristic fibrous invasions that can create conduction discontinuties and promote local circuit formation [62][63][64] . Other important differences are electrophysiologic. PV myocytes have the shortest APs of any myocyte subtype in the human heart 65 , and exhibit a clear propensity for triggered activity 66,67 . Some have suggested that a specific population of myogenic or spontaneously active myocytes exist in the PV sleeves 68 , although this has not been broadly confirmed.
A range of computational studies have incorporated some of these characteristics and begun to create a biophysically defined heirarchy of arrhythmogenic mechanisms in the PV sleeves. Earlier work from Cherry et al. 69 clearly demonstrated an important interaction of heterogeneity in myocyte-to-myocyte coupling with tissue size for promoting re-entry in a pseudo-PV sleeve. More recent computational efforts have sought to apply clinical imaging data to understand how macroscopic properties of PV electrophysiology and structure interact with www.nature.com/scientificreports/ electrical activity in the atrial free wall 12,16,[70][71][72][73][74] . These approaches often incorporate measured variations in fibre direction, together with electrophysiologic heterogeneities in regions defined at the centimeter scale, in geometries of individual patient atria 12,16,72,73 . Some have also included heuristically defined fibrotic regions within the PV sleeve 12 . Together, these studies have suggested a range of structural and dynamic substrate components are important for permitting macroscopic arrhythmia to initiate (or anchor) around the PV-LA junction, and then propagate into the atrial myocardium. However, because all of these studies employed the bidomain or monodomain models they were not able to address electrophysiologic heterogeneities or structural characteristics smaller than a length scale of several millimetres. This precludes investigating a range of properties that rely on Figure 8. EMI-based multicellular simulations of electrical propagation and re-entry in the cylinder segment of the PV sleeve tissue. We use the simulation protocol described in the "Re-entry" section for investigating re-entry. Data are illustrated for WT, and the N588K, A545P, and E299V mutations (columns). Specifically, we show the intracellular potential of the cells in the PV sleeve at five different times after stimulation (rows). Note that to improve the visual resolution, the transversal axis of the myocytes has again been stretched by a factor of 20. www.nature.com/scientificreports/ local heterogeneities or microscopic fibrotic barriers to create micro-reentrant circuits, or spontaneous foci at the PV-LA border. These proarrhythmic events have been suggested by many studies, and the PV-LA border is among the most probable locations for these proarrhythmic mechanisms in humans [62][63][64] .
Simulated pro-arrhythmic properties of AF mutations at the PV-LA junction. While the mutations we have studied here were not specifically chosen for their link to PV sleeve arrhythmias, all are channelopathies that have been causally linked to lone AF or AF that is only comorbid with other arrhythmia disorders. As such they are electrophysiologically well-characterized and known to predispose to AF without a requirement for broader confounding comorbidities (e.g. heart failure). Our implementations of these mutations faithfully recapitulated their electrophysiologic phenotypes at the level of single myocytes. Some of these were quite subtle (A545P, I to ), while others were much more severe (E299V, I K1 ). Our EMI tissue framework was able to distinguish an arrhythmic propensity in 4 of the 6 mutants, and it was only the the E375X ( I Kur ) and A545P www.nature.com/scientificreports/ ( I to ) mutants that did not exhibit clear arrhythmic activity. In all cases this arrhythmic activity was re-entrant, that is, we did not observe spontaneous activity in any of the tissue simulations. However, this likely reflects greater sensitivity of our re-entry protocol, compared to our excitability protocol, for distinguishing arrhythmic potential. In particular, we did not incorporate any instances of the PV or LA myocyte phenotypes that exhibited spontaneous depolarizations or spontaneous APs in the PV-LA tissue. Extending the variability of the myocyte phenotype to permit some instances that are spontaneously active (as would be expected of typical myocyte heterogeneity in tissue) is one direction for the framework that could be of obvious value. We did still observe marked changes in tissue electrotonic properties with our excitability simulations. These suggest that the loss-offunction c.932delC (Cx43) mutant is likely to be particularly predisposed to spontaneous focal activity, whereas I Na loss-of-function accompanying A130V may be sufficient to elicit conduction block in certain conditions. Together, these results reveal that even this initial implementation of the cell-based EMI framework is capable of assessing arrhythmic PV phenotypes involving re-entry, and that more detailed approaches to assessing spontaneous focal activity are likely possible.

Limitations.
Our results suggest that the high-resolution modeling approach permitted by EMI may improve our ability to assess the mechanisms of PV arrhythmogenesis. However, there are clearly aspects that can be improved. First, in this study we have not investigated the potential for microfibrotic regions, which are known to occur at the PV-LA border 64 , to impact conduction and introduce micro-reentry. This would require significantly larger tissue geometries, which while possible are computationally costly, and thus beyond the initial design we have applied here. Second, while we did implement random heterogeneity that was constrained by the mean properties measured in canine LA and PV myocytes 22 , it is likely that real electrophysiologic heterogeneity exceeds these constraints. As a cell-based framework, EMI is uniquely suited to applying more realistic heterogeneity, perhaps by invoking populations of myocyte models constructed from the measured variability in those same canine experiments. Third, we have not considered perturbations to calcium homeostasis in these simulations. Such perturbations are well known to play a role in PV-driven arrhythmia 75 . Fourth, we have assumed that it is sufficient to consider one layer of myocytes. The myocardial thickness of the sleeve of the pulmonary vein varies with the distance from the venoatrial junction. In 62 the thickness is reported to vary between 0.1 and 2 mm, and from 0.3 to 0.8 mm in 76 . Therefore, the use of only one layer of myocytes is an approximation selected and used in this study for computational tractability. Finally, we also acknowledge that computational cost and the relatively small computational domain used in these simulations have likely limited our ability to address important determinants of PV arrhythmia. In particular, gradients in APD from the PV ostia to the distal PV sleeve have been shown to be important for AF inducibility in this region 12 . While we have implemented local heterogeneities, these systematic gradients exist over larger spatial dimensions than our simulated geometry can allow. We also highlight two important opportunities to extend this modeling framework: (1) It would be desirable to incorporate diffusive ionic mass conservation in both the extracellular and intracellular spaces to assess longer-term changes to ionic homeostasis in these two domains-currently, the extracellular domain is assumed to be an infinite pool, thus concentrations are constant in space and time.
Significant changes in extracellular ionic homeostasis are thought to result from the high frequency AP firing that is characteristic of AF, and explicit implementation of the intracellular and extracellular spaces in EMI offers a unique basis for assessing these effects. (2) It is well known that interplay of parasympathetic and sympathetic inputs are key for inducing PV arrhythmia 66 . However, the degree to which regionalization of these inputs within the PV sleeve promotes arrhythmogenic outcomes is not well understood. Similarly, small regionalized fibrotic borders have been strongly implicated in PV arrhythmia 64 . The high spatial resolution offered by EMI could readily be applied to reveal and study the details of the arrhythmogenic potential both of these very localized effects. Table 6. Summary of the results of the EMI model simulations of the PV sleeve for WT and six mutations linked to AF. Definitions of the properties are given in the "Definition of excitability and re-entry" and "Defining biomarkers" sections, and the numbers in parentheses report the deviation from WT. The distance around the PV sleeve cylinder (of radius 1.5 cm) is approximately 4.7 cm.

Conclusions
To our knowledge, this is the first study that has simulated the PV-LA border at a resolution needed to capture mechanisms that operate at a spatial scale of a single myocyte. Our results show that the cell-based EMI framework can discriminate the functional arrhythmogenic influences of known AF mutations in simulations of the PV-LA border. We anticipate that a similar modeling strategy can be applied to investigate the role of defined local microfibrous structures, realistic heterogeneity in myocyte orientation and coupling, and clinically relevant improvement of ablation procedures in this critical PV-LA sleeve region. As such, this modeling framework provides a strong and novel foundation for assessing the unique arrhythmogenic properties of the PV sleeve.